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Numerical simulations of two-dimensional granular flows under uniform shear and external body torque were 
performed in order to extract the constitutive equations for the system. The outcome of the numerical simu- 
lations is analyzed on the basis of the micropolar fluid model. Uniform mean shear field and mean spin field, 
which is not subordinate to the vorticity field, are realized in the simulations. The estimates of stresses based 
on kinetic theory by Lun [Lun, J. Fluid Mech., 1991, 233, 539] are in good agreement with the simulation 
results for a low area fraction u = 0.1 but the agreement becomes weaker as the area fraction gets higher. 
However, the estimates in the kinetic theory can be fitted to the simulation results up to v = 0.7 by renormaliz- 
ing the coefficient of roughness. For a relatively dense granular flow (y = 0.8), the simulation results are also 
compared with Kanatani's theory [Kanatani, Int. J. Eng. Sci., 1979, 17, 419]. It is found that the dissipation 
function and its decomposition into the constitutive equations in Kanatani's theory are not consistent with the 
simulation results. 
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1. Introduction 

Collective motions of granular materials behave like fluid motions under appropriate conditions. 
However, unlike the Newtonian fluids, the basic equations for the collective motions of granular 
materials have not been well established yet. One of the difficulties of the problem lies in the fact 
that the scale of macroscopic collective motions is not well separated from the microscopic scale of 
the system such as the radius of the granular particles. Thus, applicability of arguments based on 
the scale separation would be limited. Many detailed properties of the individual particles would 
directly affect the behavior of the macroscopic flow. 

One possible way to treat such granular flows is to model them as flows of a micropolar fluid, 
a fluid with polar micro-structures such as spin 0, By applying micropolar fluid mechanics to 
granular flows, the spin of the granular particles can be coupled to the dynamics of the macroscopic 
collective motions of the granular particles. The microscopic properties of the granular particles are 
reflected in the equations of motion for the macroscopic fields through the constitutive equations, 
i.e., the relations between strains and stresses (see section [2]). 

For sparse and rapid granular flows, the equations of motion as a micropolar fluid can be 
derived within the framework of a kinetic theory. Firstly, the kinetic theory was developed without 
introducing the frictional interactions between particles and the degrees of freedom for spin by 
Savage and Jeffrey Q, and Jenkins and Savage Although Jenkins and Richman [jj, Jenkins 
and Zhang [6| and Yoon and Jenkins [3] introduced frictional interaction between particles to the 
kinetic theory, they eliminated the macroscopic degrees of freedom of the spin field by assuming 
that the macroscopic spin field is subordinate to the vorticity field. Such an assumption may be 
justified, for example, for steady flows far from the boundary. In the kinetic theories, the effect 
of frictional interactions can be absorbed into the renormalized restitution coefficient. Saitoh and 
Hayakawa [8] performed numerical simulations of two-dimensional granular flow under a plane shear 
and confirmed that the hydrodynamic equations derived from the kinetic theories agree with the 
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simulation results. In some cases, such as flows near boundaries, discrepancy between the spin field 
and the vorticity field is not negligible. The kinetic theories retaining the spin field as independent 
macroscopic degrees of freedom were developed by Lun and Savage @, Lun and Goldshtein 
and Shapiro Mitarai et al. [12| performed numerical simulation of a collisional granular flow 
on a slope and showed that the velocity and spin field profiles are in agreement with the micropolar 
fluid equations based on constitutive equations which are consistent with that in [10] . 

When the granular particles become dense enough and the volume fraction exceeds the critical 
value, the collective motions of particles stop to behave like a fluid in a sense that a finite shear 
stress is required to create an infinitesimal strain. Such a phase is called the jammed phase. The 
phase transition between unjammed phase and jammed phase is called the jamming transition. 
Scaling laws near the critical point of the jamming transition have been sug gested and verified 
in the numerical simulations by Hatano [13j, and Otsuki and Hayakawa \l4. Il5j|. The frictional 
interactions among particles were not considered in their studies. Recently, a number of results 
on the jamming transition based on numerical simulations including the frictional interactions 
have been reported (e.g., Silbert et al. [l6j], Zhang and Makse [ljj and Shundyak, Hecke and 
Saarloos [Hj|.) 

There can be a substantial intermediate regime of the volume fraction between the kinetic re- 
gion with low volume fraction and the critical region near the jamming transition. In this regime, 
interaction of n-particles with n > 2 would become important. Kanatani [l9| developed a micropo- 
lar fluid theory for relatively dense granular flows in which particles are almost regularly in contact 
with the other particles. The regime where Kanatani's theory is applicable is possibly located in 
this intermediate regime. Kano et al. [20] showed that numerical simulation of a granular flow on an 
inclined trough is in qualitative agreement with the micropolar fluid equation based on Kanatani's 
theory regarding the velocity profile. 

In this paper, we focus on the constitutive equations for granular flows. As an intrinsic nature 
of the granular flows, the spin field associated with granular particles is not subordinate to the 
velocity field of their mean flow. This situation is analogous to the case of micropolar fluids in 
which the spin field is regarded as an independent degree of freedom. As we will see in section [2j 
both the difference between the vorticity and spin fields, which will be denoted by Rji, and the 
spatial derivative of the spin field, which will be denoted by Qkjii contribute to the constitutive 
equations. From a theoretical point of view, it is desirable to analyze them separately. Therefore, 
let us consider the case of Rji ^ and Vt^ji — 0. Note that Rji ^ near the boundary. When 
sufficient numbers of particles are contained in a region near the boundary with the length scale 
smaller than the typical length scale in which the shear and spin fields changes, we may consider 
that uniform shear and spin fields (i.e. flkji — 0) with Rji ^ are approximately realized in 
the region. The situation Rji ^ and Qkji = would be also obtained by applying external 
torque to each particle. Note that, provided that the micropolar fluid picture is appropriate for the 
granular flows, the constitutive equations depend solely on velocity, spin fields and their spatial 
derivatives at the local point under consideration and independent of driving forces that generate 
the fields. A possible way to apply the external torque to each particle in experiments is to embed 
the source of angular momentum inside each particle. That is, the particle is supposed to be a 
kind of micro-machine composed of an outer shell and an inner sphere with the friction between 
them being small. Initially, the inner sphere is made to rotate with a high angular velocity by 
some means while the outer shell is not rotating. Then, the angular momentum of the inner sphere 
is continuously supplied to the outer shell through the friction until the inner sphere loses its 
substantial angular momentum. By virtue of small friction, one can realize a longer period for the 
experiment. By considering the inner sphere as an exterior system, the situation implies that the 
external torque is continuously applied to the particle (the outer shell). The inner sphere can be 
replaced by a liquid with low viscosity such as a super fluid. The actual setting of the above system 
for the experiment may be quite difficult. However, in numerical experiments, it is quite easy to 
apply the external torque to each particle. 

Taking the above into consideration, we performed numerical simulations of two-dimensional 
granular flows under uniform shear and uniform external torque field. By virtue of the external 
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torque field and the applied boundary conditions, macroscopically uniform vorticity and spin fields 
are realized and their magnitudes are controlled independently, which means that Qkji — and the 
magnitude of Rji can be controlled (see section [5]). Thus, we concentrate on the Rji dependence of 
the constitutive equations with £lkji fixed to 0. The study of the Qkji dependence of the constitutive 
equations will be the next step and will not be referred to in this paper. Unlike the preceding 
numerical studies such as [20( and [l2j, we are able to obtain not only the velocity and spin field 
profiles, which are the results of the constitutive equations, but also the constitutive equations 
directly. Since the subject of this paper is the micropolar fluid aspect of the granular flows, the 
value of area fraction is varied within the unjammed region. We compare the results from the 
numerical simulations with those from the kinetic theory by Lun (lp| , which is capable of treating 
cases that the spin field is not subordinate to the vorticity field. For the intermediate regime noted 
above, we also compared the simulation results with Kanatani's theory [19]. 

This paper is organized as follows. In section [3J a brief review of the micropolar fluid theory 
is given. In sections [3] and 2] the kinetic theory by Lun and Kanatani's theory are reviewed, 
respectively. In section [SJ comments on the two theories are given. In section [5] the results of the 
numerical simulations are shown and they are compared with the theories. In section [7] discussion 
is presented. 



2. Equations for micropolar fluid 

In this paper, we consider collective motions of particles. For simplicity, we assume that the 
particles are of the same mass m and the same moment of inertia I. Let c,(f), Wji(t) and rj(i) 
be, respectively, the velocity, the spin and the position of the particle at time t where i and j are 
coordinate indices of d-dimensional space. Here, d can formally be an arbitrary positive integer 
larger than 1. In this paper, we use the convention that the spin is expressed by a skew-symmetric 
tensor Wji whose (j, i)-th component gives the angular velocity in the (j, i) coordinate plane. Let 
F^ N \c^\ w^ 1 ', rW; • • • ; c^ N \ w^ N ' , r^; t) be the probability density function in the phase space 
of iV-particles system satisfying Liouville equation. Here, the bold letters c, w and r denote vector 
or tensor, the superscript (a) on c^ a \w^ a ' and is the index of the particle and is 
symmetrized with respect to interchanges of the particles. The s-particles set distribution function 
f( s \s < N) is given by 

/«(c« r (1 V-- ;cW,to«,rW;*)= 7 -^ TT ( [ dc^ [ dw& [ dr^ 

xF( N \ C W, W W,rW;--- ;c W t»W rW;t), (2-1) 
and the number density of the s-particles sets is given by 

n (s)( r M... ,r«;t) = f[ (J dc( Q ) J dw^ 

a— 1 * 

x /M(c«W (1) ;--- ;c««« r«;i). (2.2) 

Macroscopic fields such as the mass density field p(r, t), the moment of inertia density field pi(r, t), 
the velocity field v(r,t) and the spin field u>(r, t) are introduced as 

p(r,t) := mn (1) (r,t), Pi{r,t) := In W (r,t), (2.3) 

v(r,t) := {c) P ,t, u(r,t):= (v>) r ,t, (2.4) 

where 

(ip(c,w)) rjt := n{1) ^ r ^ J dc J dii;?/;(c,ii;)/ (1) (c,u;,r;i), (2.5) 

for an arbitrary function tp of c and w. Hereafter, indices or subscripts of spatial or time coordinates 
will be suppressed unless we need to emphasize them. 
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These macroscopic fields satisfy the following equations, 

^+^ = 0, (2.6) 

!^L +pj a t v i =Q, (2.7) 
Dv 

p-^- = dja jt + phi , (2.8) 

Pl-fyT = 2cr [j'»] + 9 k^kji + PlTji , (2.9) 

where D/Dt = d/dt+Vidi, <7ji is the stress tensor, Xkji the couple-stress tensor, bi the external body 
force and Tji the external body torque. Here and hereafter, summations are taken over repeated 
subscripts and we employ the notations, 

= T Ui) ~ ^TkkSji , (2-11) 

for an arbitrary tensor Tji. Equations (|2.6p - (|2.9p correspond, respectively, to the conservation laws 
of mass, moment of inertia, momentum and angular momentum. Note that Vi and U)ji are mutually 
independent degrees of freedom. When there is no spin field u>ji, (|2.6p - (12.8p reduce to the equations 
of ordinary fluids. A fluid with the degree of freedom of the spin field is called micropolar fluid. 
The equations of motions will be closed when the stresses Uji and Xkji are known in terms of djVi 
and ilkji '■= dkWji- The equations which yield such relations are called constitutive equations. 

Let C and W be the fluctuations of velocity c and spin w of individual particles around the 
macroscopic fields, i.e., 

C = c-v, W = w-u>. (2.12) 

The "internal energies" per unit mass e t and e r associated with the fluctuations C and W , respec- 
tively, are given by 

et:=i(aa), ^-^^(WjiWji). (2.13) 

Here, the subscripts t and r indicate the translational and rotational motions, respectively. Two 
kinds of "temperature", T t and T r , are introduced by the relations et = (d/2)T t and e r = (d(d — 
l)/4)T r . The total "internal energy" is given by eu = et + e r . Here, we have assumed that the 
particles are rigid bodies and that there is no potential force acting among particles, that is, there 
is neither contribution of elastic energy nor of potential energy to the "internal energy". One can 
show from (|2.6p - (|2.9p that the kinetic energy of macroscopic fields per unit mass, 

£k = \viVi + Y-ojjiUjji , (2-14) 



satisfies 



with 



2 4p 



p^ = y+pd i v i -$, (2.15) 



* := d k ((JkiVi) + ^<9fc (XkjiUji) + pbiVi + ^piT jl u) jl , (2-16) 

: = a {ji}Eji + &{ji]Rji + -A* ,, , (2-17) 
where we have introduced notations, 

Eji := d {j v i}l R 3 i := dyv^ - u 3i , (2-18) 
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Figure 1. Configuration of two contacting granular particles. 



and p — —aa/d is the pressure. Note that W is the work done on a fluid element per unit volume 
per unit time by the stress Oji , the couple stress Xkji , the external body force bi and the external 
body torque Tj\. From (|2.15p and the energy budget equation. 

p-^(e K + eu) =y + q-dihi, (2.19) 

we obtain 

P^- = $ - pdm + q - d t h t , (2.20) 

where q is input of the "internal energy" per unit volume per unit time, hj is the "internal energy" 
flux. The quantity $ — pdiVi gives the energy transferred from the kinetic energy to the "internal 
energy" per unit volume per unit time. When macroscopic fields p, djVi, ejj and hj are constant in 
time and space, we have diVi = from (|2.6p . and (|2.20p reduces to 

<P = -q. (2.21) 

The equation (|2.21[) implies that the energy transfer rate $ from the kinetic energy to the "internal 
energy" balances with the dissipation rate of "internal energy" — q for macroscopically steady states. 
$ is called the dissipation function since it gives the energy going out from the kinetic energy ex 
per unit time per unit volume in the macroscopically steady state. 

3. Kinetic theory for collisional granular flow 

A Kinetic theory for collisional granular flow is developed by Savage and Jeffrey Q , and Jenkins 
and Savage [4]. It is extended to include the effects of surface friction and inertial moment of 
particles by Lun and Savage and Lun From assumptions on the collisional process of two 
particles, the constitutive equations for the granular material as a micropolar fluid are derived. We 
briefly review the theory by Lun in the following. 

The dimension of the configuration space is d and the particles are assumed to be <i-dimensional 
spheres with the same radius a. The mass and the moment of inertia are, respectively denoted by 
m and I as in section [5] Let the colliding two particles be labeled as 1 and 2, and a unit vector k 
be defined by 

r (2) _ r (i) 

k := | r (2)_ r (i)| ' ( 3 - 1 ) 

(see figure Q}. In this section, we employ the notation that the quantities without (with) check " 
denote those just before (after) the contact. Let J be the impulse of the force exerted on the 
particle 2 by the particle 1 through the contact point. We have 

-(1) (1) T -(2) (2) I T fo r,\ 

Iwf = Iv$ - a{k 3 J t - kiJj), Iwf = Iwf) - aikjJ, - hJ 3 ). (3.3) 
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Let us assume that the velocity difference £j between the two surfaces at the contact point, 

6 := (c- + ak jW +){5u - kh), (3.4) 

where 

c± := c W±cW, w ± :=w^±^ 2 \ (3.5) 

changes as, 

(% - feifcjOlj = -P(6ij - kkj)£j , (3.7) 

after the collision. Here, the coefficient of restitution e and the coefficient of roughness /3 are 
assumed to be constants satisfying < e ^ 1 and —1 < (3 < 1. From (|3.2p . (|3.3p . (|3.6p . (|3.7p . we 
obtain 

_ C W — _^ 2 ) -l ^ 2 ) — _ r«- _ «„m..,»-i.. _ r^«fc.»„+ 



c- + q = - - ^kjCj ki - ^akjwt , (3.8) 

(^i c i~ ~ ^ c 7 + a(kikjW^ - kikiW+j) , (3.9) 



-(1) (1) -(2) (2) V2 

w ■■ — W ■ — W-- — w ■■ = 

3* Ji 3* 3* J( a 



where r/i = (1 + e)/2, r? 2 = (1 + /3)K/2(K + 1) and X = //ma 2 . 

Assuming the binary collisions, the equation of motion for (ip), where %j) is an arbitrary function 
of c and w, can be written as 



dt 
where 



0i(tp;r,t) -=~\J dc(1) / dc(2) / d ™ (1) / d ™ (2) / dfe 



fcc->0 



x fc-c-(2a) d fc l (^< 1 ) - V (1) ) 

x /( 2 ) (c (1) ,u; (1) ,r-afc;c (2) ,u; (2) ,r + afc;i) , (3.11) 
X{i>\r,t) :=~ y dc (1) / dc (2) / d™ (1) / d™ (2) / dfc 

fcc->0 

x fe • c-(2a) d -\^ - + ^ - V (2) ) 

x / (2) (c (1 \™ (1 \r-afc;c (2) ,™ (2) ,r + afc;i) . (3.12) 

By substituting ip = m and tp = mci in (|3.10p . one finds that the stress tensor Uji can be written 

as 

o*=o-« +<-$», (3.13) 



with 



^^-^QQ), (3.14) 



(c) ._ /) 



(mCi), (3.15) 



where <r^ denotes the kinetic contribution to the stress (jji due to the particles that cross the 

(c) 

plane perpendicular to j-axis, and a ^ denotes the collisional contribution due to the collisions of 
the two particles in different sides of the plane. 

The particle-pair distribution function is assumed to be approximated by the form 

/ (2) (c (1) , w (1) , r - ak, c (2) , w (2) , r + ak; t) 

~ 3o(2a; r, i)/ (1) (c (1) , r - afc; t)/ (1) (c (2) , w (2 \ r + ak; t), (3.16) 
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where go(r';r,t) is a radial distribution given by 

n (2) (r - ye,r + Ce;t) 



g (r';r,t) := 



(3.17) 



with an assumption that it is insensitive to the direction of a unit vector e. Assuming that a is 
sufficiently smaller than the typical spatial scale in which the amplitude of varies, we have 



(3.18) 



(3.19) 



(c, w, r ± afc; t) ~ (c, u>; r, t) 1 ± afe ■ V In (c, w. r; i) 
Let /W be written as 

/« (c, r; i) = / (1) (c, «>, r; t) [1 + <f>{c, w, r; *)] , 

where /q is the distribution function at a local equilibrium state and 4> represents the perturbation. 
The function /q is given by 

fo(c,w,r;t) 



(2wT t {r,t))i (2ivmT r (r,t)/I) 



x exp 



Ci(r,t)Ci{r,t) IW 3l {r,t)W 3l {r,t) 



2T t (r,t) 



4mT r (r,i) 



(3.20) 



Note that different temperatures T t and T r are assigned to transitional and rotational degrees of 
freedom, respectively. 

Assuming that the mean is much smaller than the typical magnitude of the fluctuation for the 
distribution of Wji, i.e., \ <C \J (m/I)T T , we have 



exp 



IW JZ {r,t)W 3l (r,t) 
4mT T (r,t) 



IwjiUji(r,t) 



2mT T (r,t) 



exp 



IWjjWjj 

'4mT x (r,t) 



(3.21) 



It is assumed that the function (f> is approximated by a linear function of degrees of nonequilibrium 
such as the symmetric and traceless kinetic stress tensor cffii , and the translational and rotational 

kinetic energy fluxes h^} and h x } given by 



PI 



{GiWijWi 



(3.22) 



In what follows, we restrict ourselves to the case of macroscopic fields djVi,ujji,Tt and T r being 
constant in space and time. In such a case, there are no energy fluxes h$ and h^- and the 
perturbation (f) depends solely on cS\ ■ Let us assume the form 



.00 



2pT t W> 1 1 ' 



(3.23) 



where the factor —l/2pT t is determined from the consistency condition that (|3.14l) is satisfied. 
Note that the total er^O is given by 



where there is no anti-symmetric part cr| k ] in the definition (|3.14j) . By substituting ip — mCiC'j/2 



-pT t Sji + a 



00 

m 



(3.24) 



into (13.10p . we obtain 



mCiCj 



(3.25) 
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where and are given by (|3.14|l and (|3.15p . respectively, and x(') by (|3.12[) . <j( c ' and 



x{mCiGj /2) reduce to functions of d 



and <j( k \ after performing the integrations with 
respect to fc,c^ and w^(i = 1,2) in p. lip and (|3.12p . Substitute the functions into (|3.25p and 
assume that the magnitudes of djVi and ujji are so small that the second or higher order terms in 
djVi and ujji can be neglected. Then we arrive at the equations, 



.0) 



2 d 



cr 



(c) 

d + 2 



1/2 



(4?/i + rf77 2 )/oa ( — ) + (2?7i + dri 2 )n k 



2 d g vpa — 772 i? 



1/2 



(3.26) 
(3.27) 

(3.28) 



where 



2( 7 rT t ) 1 /2 )0 a {2~ d - l (d + 2) + ^[677? - % + 2d mm - (d - 2)772 - 2ry 2 2 ] - |grZ/ 5 o} 



^50 



{2(3 + <%! + (d - l)(d + 3)?7 2 - 677? - 4^772 - (d 2 - 3)77! + 



(3.29) 

and v is the volume fraction which may be related to the number density n of the particle as 

n d/2 



' d+2 



-a d n. 



(3.30) 



Now, the constitutive equation for the total stress tensor a is given by 

a , ■[!■', ,■ Rji, Cl k ji = 0) = - pT t 6 ji + er[^ } + cr^-j + afy 



(3.31) 



with (|3.26p - (|3.28p . From symmetry, we have Hkji{Eji, Rji, flkji = 0) = 00- The obtained constitu- 
tive equations are a restricted version of those in [10] in the sense that constants djVi, u>ji, T t and 
T r are assumed, and they are a generalized version in the sense that d is arbitrary whereas d = 3 
in pjj. 



4. Kanatani's theory 

A micropolar fluid theory for relatively dense granular flows was developed by Kanatani (l9| . 
The theory can be outlined as follows. As in section [3J a system of the particles with the same 
radius a, mass m and moment of inertia / is considered. Let the two contacting spherical particles 
be labeled as 1 and 2. The unit vector k is the same as (|3.1I) (see figure [IJ. From (I2.12p . the 
quantities associated with particles, cf^ and (a — 1,2), can be related to the macroscopic 
fields as 

= Bi (rW) + C\ a \ wf = Uji (rW) + W$\ (4.1) 

where c\ a ^ and Wj" denote the fluctuations around the macroscopic fields. Since \r W — r ( 2 ) | = 2a 
is small compared to the typical scale in which the amplitudes of fields vary, Vi(r^) and ojji(r^) 
can be approximated by its Taylor series about x^ up to the first-order terms, i.e., 

Wi(r (2) ) = 7j 4 + 2akjdjVi , u^j(r (2) ) = Uji + 2afeAjj , (4.2) 

where we have omitted writing the argument in Vi(x^),u)ji(x^), djVi(x^) and flkjiix^). 
From (|3.4p . (|4.ip and (|4.2p . the velocity difference ^ between the two surfaces at the contacting 

1 Expand fikji{Eji, Rji, &kji = 0) in the power series of Eji and Rji. Note that the coefficient tensors are of odd 
degree. Since there is no special direction in the granular system, the coefficient tensor must be isotropic. Isotropic 
tensors of odd degree are 0. 



13401-8 



Constitutive equations for granular flow 



point is now written in terms of djVi, ujji, fijyj, ki, C± and Wj"\ Assume that the fluctuations 
Cf 3 and Wj? can be neglected and that hi is isotropically distributed. Then, we have 

/o3\ 1/2 

e:=((^)) 1/2 = (j) au, (4.3) 



where 



d „ „ 1 



2(d+2) Jl Jl 2 Jl JL 2(d + 2) 
and ((•••)) denote the average over fcj. Here, we have used 



1/2 



(4.4) 



((h)) = 0, ({kjh)) = ^ , ((hhjh)) = 0, (4.5) 



((kik m kjh)) = ^ (SimSji + SijS mi + SiiS m j), (4.6) 

which can be derived using the general expressions of isotropic tensors and kiki = 1. Suppose that 
the average energy dissipation per unit time at a contacting point is estimated by (ift;, where // is 
the kinetic friction coefficient and / is the average amplitude of the force applied in the direction k 
at the contacting point. Let iV c be a number of contacting points per a particle. Then, the energy 
dissipation per unit volume per unit time — q is given by 

2 m 

where the factor 1/2 is introduced to cancel the double counting of the contacting points. Note 
that the pressure p c due to the contacting force can be estimated by 

Pc = — — , (4.8) 

where v is the volume fraction and S is the surface area of the particle. From (|4.3I) . (14. 7[) and (|4.B|> . 
one arrives at 

-q=(2d)^ 2 n Pc 6j. (4.9) 

The estimate of p c is given as a homogeneous function of Eji, Rji and fijyj, i-e., 

p c (aEji,aRji,attkji) = a 1 * p c (Eji,Rji,Q).ji) with a constant The form of the function depends 
on whether the flow is slow or fast, and here we omit the review. See the original paper \m for 
the details. 

Let us restrict ourselves to macroscopically uniform and steady states. From (|2.2ip and (|4.9p . 
one obtains 

$ = {2d) 1 ^ t ip c uj, (4.10) 

which is a homogeneous function of Eji, Rji and fijyi, i-e., <fr(aEji, aRji, a£lkji) = 
a^&(Eji, Rji, ilkji) with £ = £' + 1. With the help of Euler's homogeneous function theorem, 

•-cl^+^+ass *)' (4 - n) 

the choice 

1 » 1 » 119$ 

~CdE~i' ° m ~ C ' 2 Mfc ^ ~ ' ( ] 

is made for constitutive equations which are consistent with (|2.17l) . 
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5. Comments on the kinetic theory and Kanatani's theory 

In the kinetic theory by Lun, the binary collision is assumed, i.e., n-particle collisions for 
n > 2 are neglected, and the particle-pair distribution /( 2 ' is assumed to be approximated by a 
product of go and /W as (|3.I6p . These assumptions seem to be valid for small volume fraction, 
v <C 1. This is because, in such a case, n-particle collisions for n > 2 would be rare and the 
velocities of the two particles before a binary collision would be statistically almost independent. 
Furthermore, macroscopic fields djVi and u)ji are assumed to be small in comparison to their 
microscopic fluctuations. All of these assumptions would become inappropriate with the increase 
of v. 

Even when v is small, there are a few points one should consider carefully in the kinetic theory 
by Lun. Regarding the distribution function /q in (|3.20p . different temperatures T t and T r are 
assigned to transitional and rotational degrees of freedom. This is not a redundant setting since the 
violation of equipartition of energy between transitional and rotational degrees of freedom indeed 
occurs. For example, Huthmann and Zippelius [2~D | showed that the equipartition is immedi ately 
destroyed even if it is initially satisfied. Some studies (e.g., Goldhirsh, Noskowicz and Bar- Lev [22j, 
Kranz et al. [23]) suggest that there are substantial deviations from Gaussian and correlation 
between C, and Wji in the distribution function /W. Therefore, there can be a substantial deviation 
from (j3~Ti?|) with fX2H)l and in /W. Since our aim is to extract information for the constitutive 

equations, full detailed information on /W is not necessary. However, we do not know how much 
information is sufficient for our purpose a priori. Here, it is assumed that the effect of the corrections 
on moments higher than two in /q is not so large. 

Note that it is suggested by a number of studies that (3 depends on the angle d := arctan(£ t /£ n ') 
of the oblique collision of the particles, where £„ = ki£i and £t = y£i£i — Walton and Braun [24| 
derived 

f-l+(l+i)M(l + e)cottf (§>§ ) , v 

\/3o (<> < 0o) ' 

where j3o is the maximum value of the coefficient of roughness and i9o is the critical angle. The form 
of j3 is verified in the numerical simulation by Saitoh and Hayakawa [8]. Therefore, the constant f3 
assumed in the kinetic theory by Lun, should be regarded as an approximation. Again, we do not 
know a priori how much detailed information on /3 is required to obtain constitutive equations. It 
should be verified whether the simplified model of constant (3 is adequate herein. 

In Kanatani's theory, a relatively large v is considered. Every particle is in contact with some 
other particles during almost all the time. The estimate of the energy dissipation per unit 
time per a contact point implies that the velocity difference £ is maintained to the order of the 
macroscopic time scale. However, it is not clear whether this is the general case for a large v. 
The velocity difference £ possibly decays in the order of microscopic time scale during the contact 
due to friction. Furthermore, among many choices of constitutive equations that are consistent 
with a given dissipation function &(Eji, Rji,Vtkji), a particular choice (14. 12(1 is made. When the 
order of homogeneity £ is 2, the choice of constitutive equations (|4.f 21) leads to the linear response 
satisfying Onsager's reciprocal relations [25| . Thus, the choice is valid. However, when C 7^ 2 and 
the constitutive equations are nonlinear, such a validation is not possible. 

Some of the points given above will be examined by means of numerical simulations in the next 
section. 



6. Numerical simulations 
6.1. Setting 

We performed numerical simulations of a system of granular particles using a distinct element 
method (DEM). In DEM, the interaction forces for every contacting pair of particles are calculated 
at each time step, and the equations of the motion are solved using a difference method. The 
dimension d of the configuration space is 2. All the particles are disks of the same radius a. The mass 
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is uniformly distributed inside the particles and thus the mass m and the moment of inertia / of the 
particles are related as K = I/ma 2 = 1/2. As a mechanism of the contact, we assume Hooke's law 
of elasticity in the direction parallel to k of (|3.ip , that is, F — nl where I = max(2a— ~x^\,0) 
is the overlap length and n is the force constant. The kinetic friction force [j,F is applied at the 
contacting point in the direction to reduce the velocity difference £ between the surfaces of the two 
particles (see section 0]). We have chosen the above model of contact since it is one of the simplest 
models that induces both energy dissipation and rotational degree of freedom of particles, that 
are essential elements to consider the micropolar nature of the granular flow. For simplicity, we 
have neglected nonlinearity in the relation between I and F, inelasticity in the direction parallel 
to k and the static friction in the direction of the velocity difference £. By virtue of simplicity of 
the model, the microscopic characteristics of granular particles are completely determined by four 
parameters, the radius a, the mass m, the force constant k and the kinetic friction coefficient \i. 

N particles are put in a square domain with the length of sides L x and L y . The area fraction 
v is given by v = Nira 2 / L x L y . We fix L x /a — L y /a — 200 for all the simulations except for the 
case v = 0.1 in which L x /a — L y /a — 600. The time scale associated with the collision is given 
by i C oi = (mj n) 1 !' 2 . The time step At of DEM is set to 10 — 2 t co i, which is sufficiently smaller than 
t ca \ so that the collision process is resolved in the simulation. Since there is no damping in the 
direction parallel to k, the time scale associated with the damping is idamp = oo. 

For the kinetic friction coefficient, we used n = 0.3 and 0.8. Note that kinetic theories which 
eliminate the spin field as an independent field are suggested 0, 0] for small kinetic friction coef- 
ficient. Therefore, the aspect of granular flows as a micropolar fluid would become significant for 
large /i. In order to investigate this aspect, we have chosen somewhat larger /i compared to those 
in some measurements, e.g. /i < 0.2 in [26j |. 

We employ the Lees-Edwards periodic boundary conditions (2?| for the velocity c of particles, 

Ci(x = L,y,t) = Ci(x = 0,y,t), (6.1) 
c t (x,y — L,t) — Ci((x - Lyjt) mod L x ,y = 0,t) + LyjSix , (6.2) 

where 7 is a constant. The external torque fj% = f(Sj y Si X — 5j X Si V ), constant in time, is applied to 
every particle. 

The macroscopic fields p(x,y), Vi(x,y) and ujji(x,y) are defined by the spatial averages, 

a a 

where c^ Q ' and are, respectively, the velocity and the spin of the particle labeled by a, := 
^2 a 1 and ^2 a denotes the summation over the particles such that x — Ax/2 ^ x^ < x + Ax/2 
and y — Ay/2 ^ < y + Ay/2. The notations Eji(x, y), Rji(x, y) and fijyj(a;, y) are introduced 
similarly to (|2.18[) . The simulations were performed up to the time that the velocity field relaxes to 
a quasi-steady state. For all the sets of parameters that we investigate in this paper, the velocity 
field relaxed to a uniform simple shear profile with the shear rate 7, i.e., Vi(y) = (jy + A)8i X where 
A is a constant. Here, we have chosen Ax — L x so that vt is independent of x. In terms of Eji(y), 
we have 

E H (v) = \ {8jv S ix + SiySjx) ■ (6.4) 

It was found in the simulations that, when Eji{y) relaxed to a uniform field, u)ji(y) also relaxed to 
a uniform field. Consequently, we have 

Rji{y) = R(S jy S. lx ~ S iy S jx ) , flkji(y) = 0. (6.5) 

The shear rate 7 is directly controlled by the boundary conditions (|6.1[) and (16. 2p . We observed 
R = 0, i.e., the spin uiji is subordinate to duvs, when there is no external torque f . The value of 
R is controlled by changing the value of f . In the present study, we varied f in the range f ^ 0. In 
such cases, we observed R ^ 0, i.e., Lo yx ^ di y v x ], at the quasi-steady states. 
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Figure 2. Stresses a± normalized by mj 2 as functions of Rf-y for the mean shear 7t co i = 
0.5 x 1(T 3 , 1.0 x 10 -3 and 2.0 x 10 -3 . The friction coefficient fi and the area fraction v are fixed 
to fi = 0.3 and v — 0.7. 



Let us consider a line Xi — h where ft. is a constant. In the simulation, the kinetic contribution 
to the stress tensor, cjj^, is estimated by 

J a 

where the summation ^ s taken over the particles a which cross the line Xj = h during the 

time step At. The contribution of contacts to the stress, cr^f, is estimated by 

^W = ^£ (/ V^ (6-7) 

where is the contacting force applied to the particle /3 by the particle a and the summation 

is taken over the contacting pairs {a,/3} such that Xj ^ h < Xj- The mean stress <7ji(h) 
over a line Xj — h can be averaged over h to give the mean stress <jji of the whole domain, 



with 

r (k) 



^ m »! B) « , i w . ( 6 - 9 ) 

y a 

x-E^r"^-^)' (6 - 10) 



where the summations ^ Q and ^ Q a are, respectively, taken over the particles and the contacting 
pairs. Note that the pairs {a,/?} and {/?,cv} are not distinguished and are counted only once in 
the summation ^ Q ^. 

Under the conditions (|6.4[) and (j6.5|) . the constitutive equations read 

a± = a ± (i,R), (6.11) 

where we have introduced the notations er + := <?{ yx } an d G - '■= a [yx\ - There are three time scales 
involved in the dynamics of the present system. They are r= I7I -1 , £r := and the time 

scale of collision t co \ = (m/k) 1 / 2 . In the simulations, we set tj ~ tn ^> t co \. Let us assume that, in 
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Figure 3. (a) The value of the radial distribution function at r — 2a, go = g(2a), in the 
simulations for various v and R = 0. Dashed line indicates the value of the low density expansion 
to the first order, (b) The normalized pressure p/mj 2 and the go as functions of \8v\ — \v — 
where vj = 0.84. The slope of the lines indicate the scalings in theory in [15|. 



such a case, a+ and cr_ do not depend on t co \. Then, a dimensional consideration yields similarity 
forms. 



n (~..R) = su,n ( -. |-Vt ( = rn\y A /a :: 



(6.12) 



where cr± are dimensionless functions, which may depend on the kinetic friction coefficient /j. 
and the area fraction v. The origin of the function sgn(7) is a requirement of the symmetry 
<j±{— j,— R) = — £+(7, R). Note that, in the case of R = 0, (|6.12p reduces to the well known 
Bagnold scaling [28], a + oc j 2 . In figure [51 u±/m^ 2 from the simulations are plotted as a function 
of R/i for jt co \ = 0.5 x 10~ 3 , 1.0 x 10 -3 and 2.0 x 10~ 3 (/x and v are fixed). The collapse of the data 
implies that the similarity forms (|6.12p are valid within the parameter range of the simulations. By 
virtue of (|6.12p . we can fix 7 and vary only R to obtain the constitutive equations. In the following 
simulations, we put ^t co \ = 1.0 x 10~ 3 . 

The value of the radial distribution function at r — 2a, go = go (2a), in the simulations is given in 
figure EJa). The data are given for R — 0, but they are almost independent of R. Note that, for the 
elastic system (e = 1) without friction (f3 = —1) and shearing (7 = 0), the low density expansion of 
go at a thermal equilibrium state can be given up to the first order of v as go = 1 + ^[(8/3) — (v3/ n )] 
(see, for example, [29]), and it is also shown in a dashed line in figure[3]for comparison. One can see 
from figure |21 that go of the simulation is approximated by the low density expansion for relatively 
low area fraction v = 0.1, but it largely deviates from the low density expansion when v 0.6. It is 
known that, when v exceeds a critical value vj, which is called a jamming point, the system enters 
the jammed phase, in which the shear stress remains nonzero in the limit of zero strain, i.e., yields 
stress. When v approaches i/j, quantities such as the pressure p and go are expected to obey certain 
scaling laws. It is found that the present system suffers the phase separation between crystallized 
region and fluid region for v > 0.815, before entering the jammed phase. Since it is impossible to 
estimate i/j in the present system, we borrowed the value i/j — 0.84, which is a round-off value of 
vj reported in [lij for two-dimensional poly-disperse frictionless granular system, as a reference. 
In figure [3fb), p and go are plotted as functions of \6v\ — \v — z/j| for various v and R = 0. Since 
the crystallized phase is out of the scope of the present paper, data for v > 0.815 are omitted in 
the figure. The scaling laws p ~ |oV|~ 4 and go ~ predicted by Otsuki and Hayakawa [3] for 
the particles with frictionless, linear spring interactions, are indicated by lines in the figures as a 
reference. Since there are linear spring and frictional interactions in the present model, the scaling 
law could be modified. However, since the present system experiences crystallization, the scaling 
range is very narrow, in case it exists. 
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For each particle-pair collision, we can define the interfering number nj nt , which is the maximum 
number of other particles contacting with the particle-pair simultaneously during the collisions of 
the particle-pair. n- m t = implies that the particle-pair collision is not interfered by the other 
particles. In figure^ the probability P(rii nt ) is given for v — 0.1, 0.7 and 0.8 in the case of /j, = 0.3. 
The data for /j, = 0.8 are omitted since they almost collapse with those of fi — 0.3. The fact that 
P(rtint) has a peak at ni nt = for v = 0.1 and 0.7 is consistent with the kinetic theory. However, 
the peak shifts to n- mt — 2ioxv — 0.8. This suggests the failure of the applicability of the kinetic 
theory based on binary collisions. Note that the coordination number Z at v = 0.8 is about 0.6 — 0.7 
and it is still much smaller than the value Z = d + 1 = 3 for the isostatic condition [30j • 

Summarizing the above observations, we can consider that 0.7 ^ !/ ^ 0.8 is located in an 
intermediate regime where v is neither low enough so that the kinetic theory is applicable nor it is 
high enough so that the system is quite near the jamming point vj. 

The above conclusion might seem to be inconsistent with the results of the numerical simulations 
of frictional granular shear flows by Otsuki and Hayakawa [31]. According to when friction is 
introduced to the system, the critical point vj of the area fraction splits into two points vs and 
i^L, where the former belongs to the solid branch and the latter liquid branch. The point v = 0.80 
for n = 0.8 is located in the region v$ < v < Vi,, which implies that the selection of the solid 
or liquid phase depends on 7t C oi- Although the figures for 7i co i = 10~ 3 , which is the value in 
the present study, are not given in (3l| . it can be located in the solid (jam med) phase. Such a 
discrepancy between the results of the present simulations and those in [31| can occur since the 
details of the system are different in the present study and [3l[. For example, the radius of par ticles 
is unique in the present study whereas particles with four different radii are present in [31|; the 
normal viscosity is not introduced in the present study whereas the normal viscous constant is 
presumably non-zero (although the actual value is not given) in 3l|; and contacting particles are 
always slipping in the tangential direction in the present study whereas contacting particles can 
stick when the normal contact force is strong enough in Especially, the latter fact that the 
particles are always slipping in the present study would act to raise the value of vj (z/g or Vi) from 
those in 13 ill. 



6.2. Comparison with the kinetic theory 

In order to compare the simulation results with the kinetic theory by Lun, we need to determine 
the coefficient of roughness /3 in the simulations. In the simulations, we define fi by 

(Sij - kikj)i 3 = -0(Sij - kikj)£j , (6.13) 
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Figure 5. Averaged coefficient of roughness as a function of the oblique collision angle i? for 
various ^ and /i. The solid line indicates the theoretical function (J5TTJ. The dashed line shows 
P($) = cosi?, the form for the random collisions without spin. The dot-dashed line shows P(i?) 
of the modeled collision with ? = 0.5. 
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Figure 6. Averaged coefficient of roughness j3 obtained from the simulations (circle symbols) 
and the coefficient of roughness /3 e ff fitted to the simulation data on the basis of the kinetic 
theory by Lun (square symbols) for various v and [i. The error bars indicate A/3. 



where k is the unit vector parallel to — fW. Since a collision has a finite duration in the 
simulations, k and k are distinguished here. The particle-pair possibly contacts with other particles 
during their own collision, especially when the area fraction v is high and the duration of the 
collision is long. is defined as the average of (3 over the particle pairs with the oblique collision 
angle i?. In figure [5j /3(f?) is plotted with the probability density function P(fi) of i?. The data is 
for the simulation without the external torque and the dependency of on the external torque 
is weak (the figure omitted). For all (y,\i) we investigated, we have —1 ^ (3(d) 1 which is 
consistent with the assumption — 1 ^ (3 ^ 1 in the kinetic theory. When v = 0.1, /3 is in agreement 
with the theoretical estimate (|5.ip with /3rj = 0. As v increases, /3($) deviates from the estimate 
(|5.1[) . The deviation is in the direction of decreasing the ■& dependence of (3(d). This deviation can 
be understood as the effect of the collision of the particle-pair with the other particles during the 
collision within the particle-pair. The other particle randomly changes the velocity difference £ just 
after the particle-pair collision. Therefore, j3 becomes a random variable for fixed $. For reference, 
the standard deviations A/3 at <8 = 1 are 0.004(i/ = 0.1),0.05(i^ = 0.7),0.31(^ = 0.8) for /i = 0.3 
and 0.0008(^ = 0.1),0.05(^ = 0.7),0.14(i/ = 0.8) for /i = 0.8. As expected, A/3 increases with v. 

In figure[5j the dashed line shows cos the form for the random collisions without spin. Particles 
tend to collide with large angle d when spin and frictional interaction are introduced. This can be 
understood by a simple model of a particle-pair with the same spin w colliding with each other 
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with the relative velocity c. Let b be the impact parameter and 9(— it/2 ^ 9 ^ the angle 

between c and k. Then, we have 

& = 2asin9, tani9 = tan9 — , (6.14) 

cos 9 

where q := 2aw/c. Assuming the uniform distribution of b in —a ^ b ^ a, we have for the 
probability density P(i?) of i? {-it/ 2 ^ d ^ it/ 2) as 

P(tf) = ~ 1 , (6.15) 

2(1 + tan 2 6)(vT+ tan 2 9 - c tan 9) cos d 

with 



tan9= "-rsv-"* ( 6 .i 6 ) 



tantf + <^\J tan 2 i? + (1 - <j 2 ) 

In figure [SI the symmetrized probability density function P(i?) := P($) + P(— 1?)(0 ^ "d ^ tt/2) 
with q — 0.5 is given in a dot-dashed line for reference. The model probability density function 
roughly fits the corresponding simulation results for v = 0.1 and 0.7. In this simple model, c and 
uu are represented by constants and their probability density functions are not considered. If we 
take c and w as the root mean square of the fluctuations in the simulations, we have q ~ 0.7 and 
it is not far from <^ = 0.5. Thus, P(i3) for v ^ 0.7 can be roughly explained by the effect of spin. 
However, it is found that P("d) of v = 0.8 is hard to be fitted by that of the model. Especially, 
the bump around i9 = is not explained by the model. It seems that n-particles interaction with 
n > 2 should be taken into account in order to understand P($) for v = 0.8. 

In the kinetic theory by Lun, f3 is treated as a constant. Whether we can, or cannot, ap- 
proximate dependent /3(9) by a constant is not obvious. Here, we take a simple average, 
j3 := f* /2 dtf/3(tf)P(tf). The values of (3 for various v and fi are given in figure [5] 

The stresses cr+ , o+ , and cr_(= ci ') normalized by mj 2 as functions of P/7 obtained in the 
simulation are given for various v(— 0.1, 0.7 and 0.8) and fi(= 0.3 and 0.8) in figures [7] and [5] The 
constitutive equations (|3.26[) . (|3.27p and (|3.28p from the kinetic theory are also given in the figures. 
Here, the values obtained in the simulations are used for the "temperatures" T t and T r , and the 
value of the radial distribution function at r — 2a, go — go{2a). The coefficient of restitution e is 
set to 1 following the situation of the simulation. f3 is set to be /3. 

When the area fraction is as small as v = 0.1, the kinetic theory is in good agreement with the 
simulation results with regard to o±) where the relative deviations are within 5%. The relative 
deviations for is about 20%. However, since is much larger than \<r^\ 7 the deviations 

for a + (= erf + cr^ ' ) are almost determined by those of erf) . The agreement between the kinetic 
theory and the simulation for cr_ is not as good as that of <r+. The relative deviations are about 
40% for /J, — 0.3 and 20% for /.i = 0.8. This may be due to the approximation of the constant {3 in 
the theory. A better agreement for fi = 0.8 compared to that for /i = 0.3 can be explained by the 
fact that the range of i9 in which /?($) can be approximated by a constant is wider for fi = 0.8. 

As v increases, the ratio a^y /a+ ^ increases. We have cr+Vcf ~ 20 for v = 0.7. In spite of 
this high ratio, the kinetic theory still gives fairly good estimates of erf and er^ ' whose relative 

deviation from the simulation results are about 20% for erf and 10% for <j^\ As expected from 
the fact that the peak of P(n- m t) is shifted from nj nt = for v = 0.8 (see figure [4]), the kinetic 
theory is no more appropriate for v = 0.8 and the formal application yields the overestimates of 
by about 40% or higher. For er_, the overestimate is larger by the factor of more than 2 (the 
figures omitted). 

Now, let us find the effective /3 e ff that fits best with the constitutive relations obtained in the 
simulation. For every set of (/i, u), we chose (3 = f3 e g where (3 c g is the value of [3 which minimizes 
the summation of squared relative errors, i.e., 

f laftRil^-ai^iRi)} 2 k th (P 4 ; 6) - a sim (Ri)] 2 ) , 

c ^ - ? 1 + g^r ' (6 - 17) 
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Figure 7. Stresses cr+\ [(a),(c)], and — <r_ [(b), (d)] normalized by m-y 2 as functions of 

—R/j for the area fractions v — 0.1 and 0.7. The kinetic friction coefficient /i = 0.3 for all the 
figures. Shown in (a),(c) are and in the simulation (white circle and square), in the 
kinetic theory 1)3.261) . (|3.27p with /3 — /3 (black circle and square, interpolated with line) and in 
the kinetic theory with /3 = /3 c ff (vertical and diagonal cross). Shown in (b),(d) are <j_ in the 
simulation (white circle), in the kinetic theory p.28[) with /3 = f3 (black circle, interpolated with 
line) and in the kinetic theory with /3 = (3 e g (vertical cross). 

with the superscripts 'th' and 'sim' denoting the kinetic theory and the simulation respectively, 
and the subscript i the index of simulations with different values of R. We excluded i with Ri = 
in the summation for ct_ since ct_ (i? = 0) is ideally and thus it is not appropriate to use the 
relative error. In the vicinity of /3 = /3 e ff, the function C(/3) can be approximated by a quadratic 
function of (3, i.e, C(/3) ~ (7(/3) — a((3 — (3 c tf ) 2 + C m i n . We determined the error A/3 of to satisfy 
C(P±A0) = 2C aSn . 

The obtained values of /3 e ff and A/3 are given in figure [5] together with /3. Stresses <7+\a^ 
and tr_ obtained from (|3.26p , (|3.27p and (|3.28p with /3 = (3 c g are given in figures [7] and [SI One 
can see from the figures that ct_ can be well fitted to the simulation data by adjusting j3. On the 
other hand, ct+ is less sensitive to /3 and so the deviation from the simulation data is not effectively 
reduced by adjusting /3. The values of /3 e g are considerably smaller than /3 for all v and /j, that we 
have investigated. The discrepancy between /3 and j3 e s becomes larger as v becomes larger. 

6.3. Comparison with Kanatani's theory 

As we have confirmed in the previous subsection, the estimates from the kinetic theory are in 
fairly good agreement with the simulation results up to v — 0.7. When v = 0.8, the peak of P(n; nt ) 
is not rijnt = 0, that is, most pair collisions are interfered by the other particles. Consequently, 
the kinetic theory is not applicable. In this subsection, we examine the applicability of Kanatani's 
theory for relatively dense granular flows to the results of the simulation for v = 0.8. 
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Figure 8. Stresses 0"+', 0+^ [(a),(c)], and — cr_ [(b), (d)J normalized by mj 2 as functions of 
—R/j for the area fractions v — 0.1 and 0.7. The kinetic fraction coefficient /1 = 0.8 for all the 
figures. Shown in (a),(c) are 0+ and in the simulation (white circle and square), in the 
kinetic theory (|3. 2611 . (I3.27|l with ft = f3 (black circle and square, interpolated with line) and in 
the kinetic theory with ft — /3 c g (vertical and diagonal cross). Shown in (b),(d) are a- in the 
simulation (white circle), in the kinetic theory (|3.28|l with /3 = P (black circle, interpolated with 
line) and in the kinetic theory with ft = f5 e s (vertical cross). 

500 r 1 1 1 1 1 




(L5 1 , L5 2 2.5 

-Rft 



Figure 9. Dissipation function $ normalized by my 3, as a function of —R/j in the simulation 
for fi = 0.3 (white circle) and /i = 0.8 (white square). Dashed lines are fitted functions of the 
form (|6.20p . Also shown are the normalized dissipation function (|4.10p in Kanatani's theory with 
the value of p c from the simulation for fi — 0.3 (black circle) and fi — 0.8 (black square). As a 
guide, interpolations of the calculated data points are given in solid lines. 
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Figure 10. Normalized stresses a±/m'y 2 as functions of — R/j in the simulations and those based 
on the decomposition by Kanatani (|6.19l) for /j, — 0.3 [(a)] and /i = 0.8 [(b)]. v — 0.8 for both 
figures. Circle and square symbols denote normalized <r+ and cr_ in the simulation, respectively. 
Thick and thin dashed lines are the fitted lines H6.27[) . (|6.28[) to cr+ and cr_, respectively. Thick 
and thin solid lines indicate the decompositions (I6.25I) . (|6.26p for cr+ and <r_, respectively. 



Under the conditions (|6.4I) and (16.5I) . the dissipation function <f> of (12 . 1 7[) reduces to 

m,R) = a + (i,R)j + 2a-('y,R)R. (6.18) 

In figure O $ obtained from the simulation is plotted with the estimate (|4. 10(1 of $ in Kanatani's 
theory. In the estimate, the values from the simulation are used for p c . One can see from the figure 
that they are not in agreement. Note that the estimate (|4.10[) depends linearly on /i besides the 
possible n dependence in p c . However, <f> from the simulations depends on /i less sensitively. It 
is suggested from the disagreement of the [i dependence that the process of energy dissipation 
assumed in Kanatani's theory is not valid in the present situation of the simulation. Although, at 
v = 0.8, the effect of a multi-contact of the particle may be significant, the velocity difference at 
the contact point would not be completely sustained during the contact. 

Now, let us examine the choice (|4.12l) of the constitutive equations in Kanatani's theory. In 
order to focus on the examination of the choice (|4.12[) . let us use $ obtained in the simulation 
rather than $ of the estimates in Kanatani's theory. The form f)6 . 1 8[) with (|6.12p implies that the 
degree of homogeneity ( is 3. The choice (|4.12l) reads, in the present context, as follows: 

19$ 1<9$ .„ . 

The normalized dissipation function $ := $/m7 3 in the simulation can be fitted by a polynomial 
of r := R/j, 

$ = A + Br 2 , (6.20) 

with 

A ~ 45, 5-9.4, O = 0.3), (6.21) 
,4-52, £-15, O = 0.8), (6.22) 

(see figure [9]). Then, the choice f|6. 19[) implies that 

d+= A+ + B + r 2 , a-=A_r, (6.23) 

with 

A+ = A, B+ = ~, A - = J- (6-24) 
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From l(6T2lj) . ([6~22| and (j!T2"4l . one obtains 



~ 45, 
A+ ~ 52, 



B + ~3.1, A 
5+ ~ 5.0, A 



3.1 



5.0 



fa 
fa 



0.3), 
0.8). 



(6.25) 
(6.26) 



On the other hand, we obtain 



A+ ~ 45, 
A+ ~ 52, 



B+ ~ 1.0, A 
5+ ~ 2.3, A 



4.2, 
6.3, 



(/<■ 
(/<■ 



0.3), 
0.8), 



(6.27) 
(6.28) 



by fitting the function of the form (|6.23l) directly to the simulation data. One can see that they 
disagree. Especially, in the decomposition based on (|6.19p , we have B + = A- according to (|6.24p . 
However, this is not the case in the simulation. The disagreement can be also checked in fig- 
ure [TUJ the normalized stresses a± = a-±/m A f 2 of (|6.23p with (|6.25p and (|6.26p are plotted together 
with the normalized stresses directly measured in the numerical simulation and their fits (|6.23p 
with (|6.27p and (|6.28p . Thus, the choice (|4.12p of the constitutive equations made by Kanatani is 
not appropriate in the present situation of the simulation. 

7. Discussion 

We confirmed in section R)T2"l that the constitutive equation for er + in the kinetic theory by Lun is 
in good agreement with the simulation results for a small area fraction v = 0.1. But the agreement 
for cr_ is not as good as tr + . Since er_ is sensitive to /3, this may be due to the oversimplification of 
modeling j3 by a constant in the Lull's theory. As v increases, discrepancy between the kinetic theory 
and the simulation results on both <r + and <r_ increases. However, cr_ in the kinetic theory can be 
formally fitted to the simulation results by introducing the effective coefficient of roughness /3 e ff 
which is substantially smaller than the actual averaged beta j3 in the simulation. This implies that, 
although the particle-pairs almost get stuck {fi ~ 0) in the tangential direction during the collision, 
the macroscopic field feels as if the particles were slipping (0 e f£ ~ — 1)- A possible scenario is that 
the collective motions of some stuck particles affect the macroscopic field as motions of virtual 
particles with renormalized mass m, radius a, coefficient of restitution e, and ft. In the present 
study, we renormalized j3 with fixed to, a and e. Then (3 is renormalized to reduce its value. 

We have seen in section l6~3l that Kanatani's theory is not applicable to the present situation of 
simulation (y = 0.8) in the following sense; (i) the kinetic friction coefficient /i dependence of the 
dissipation function $ and (ii) the choice (|4.12p made for determining the constitutive equations 
from <E>, are in disagreement with those in the simulations. 

To summarize, there is a regime of relatively dense (y ~ 0.8 in the present simulation) granular 
flows in the form of micropolar fluid, where physical pictures based on neither kinetic theories 
nor Kanatani's theory are adequate. In this regime, a new microscopic description of particle 
interaction, that is different from mutually independent short-time collisions in the kinetic theory 
or long-time contacts with sustained velocity difference in Kanatani's theory, is necessary. As we 
have seen in sections l6.1l and [6T2l there are considerable changes in the interfering number rii n t, /3(#) 
and P(i?) when v increases from 0.7 to 0.8. These changes suggest the significance of the effect 
of n-particle interactions with n > 2 in this regime. It would be a future study to see whether 
the effect is included in the kinetic theory with renormalized parameters as we have preliminary 
analyzed or in some variant of Kanatani's theory with an appropriate dissipation function and its 
decomposition. 
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MaTepia/ibHi pIbhjihhji p,j\9\ rpaHy/ibOBai-ioro noTOKy 3 oflHopiflHMM 
cepeflHiM 3cyBOM i cniHOBMMM no/iflMM 

K. TaKeni, K. l/louiifla, T. ApiiviiMy 

Bnm,a LUKO/ia cpyHflawieHTa/ibHux Ta npuK^aflHi/ix HayK, ymBepci/iTeT m. L(yKy6a, l6apaKi, JlnoHia 

fl/lfl TOTO, U\0& BMOKpeMWTM MaTepiajlbHi piBHflHHn fliin CHCTeMH 6yjlM 3fli£lCHeHi HMCflOBi CMMy/lflLn'T flBOBI/1- 

MipHwx rpaHyjibOBaHwx noTOKiB nifl flieio oflHopiflHoro 3cyBy i 30BHiujHboro KpyTWJibHoro MOMeHTy. Pe3y.nbTaT 
Mwce/ibHwx cwMyjiaLiiw npoaHa/ii30BaHO Ha ocHOBi MOfle/ii MiKponojiapHoro n/iwHy. B cwMyjiau,mx pea.ni3ye- 
Tbca no/ie oflHopiflHoro cepeflHboro 3cyBy, nKe He e niflnopnflKOBaHe no/ito BwxopoBOCTi. Ouj h kh Hanpy>Keb, 
3po6jieHi Ha ocHOBi KiHeTWHHoT TeopiT JlKDHa [Lun, J. Fluid Mech. 233(1991) 539], flo6pe y3rofl>KyK)Tbca 3 
pe3yjibTaTaMW cwMyjiaujM b o6jiacTi HH3bKHX nacTOK v = 0.1, ajie y3rofl>KeHHiR noripwyeTbca, flKmp 14a Be/11/1- 
hm Ha 3pocTae. ripoTe, oujHKM, 3po6jieHi b KmeTi/iHHiH TeopiT MO>KyTb 6y™ niflirHam flo pe3y/ibTaTiB ci/iMyjinLn'ti 
a>K flo v = 0.7 w/iaxoM peHopMa/iisaajT KoecpiujHTa ujopcTKOCTi. j\s\9, BiflHOCHO rycToro rpaHyjibOBaHoro no- 
TOKy (y = 0.8), pe3yjibTaTM cwMyjiaujti Ta ko>k nopiBHioioTbCfl 3 Teopieio KaHaTaHi [Kanatani, Int. J. Eng. 
Sci 17(1979) 419]. 3HawfleHO, mo flwcwnaTWBHa cpyHKujn i TT fleKOMno3i/m,i:R b MaTepia/ibHi piBHAHHn b TeopiT' 
KaHaTaHi He y3rofl>KyeTbC!R 3 pe3yjibTaTaMi/i ci/iMynsujCi. 

K/iiOHOBi csiOBa: rpanynbOBaHiAiA noTiK, MarepiajibHi p'ibhhhhh, MiKpononxpHtAi/i nniAH, KiHertmHe piBHHHHX 
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